{
    "cells": [
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "import numpy as np\n",
                "import matplotlib.pyplot as plt\n",
                "plt.rcParams['figure.figsize'] = 10, 6\n",
                "np.random.seed(32371)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# Validation for hyperparameter selection for ridge regression\n",
                "In this notebook, we will see how to use validation to select the optimal regularization coefficient for ridge regression, working on a synthetic dataset so as to compare our observations to the theoretical optimum.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "def featurize_fourier(x, d, normalize = False):\n",
                "    assert (d-1) % 2 == 0, \"d must be odd\"\n",
                "    max_r = int((d-1)/2)\n",
                "    n = len(x)\n",
                "    A = np.zeros((n, d))\n",
                "    A[:,0] = 1\n",
                "    for d_ in range(1,max_r+1):\n",
                "        A[:,2*(d_-1)+1] =  np.sin(d_*x*np.pi)\n",
                "        A[:,2*(d_-1)+2] =  np.cos(d_*x*np.pi)\n",
                "\n",
                "    if normalize:\n",
                "        A[:,0] *= (1/np.sqrt(2))\n",
                "        A *= np.sqrt(2)\n",
                "    return A\n",
                "\n",
                "def generate_X(n, d):\n",
                "    \"\"\"\n",
                "    This function generates a random n x d synthetic data matrix X such that X^T X = I\n",
                "    \"\"\"\n",
                "    x = (np.random.uniform(-1, 1, n))\n",
                "    return featurize_fourier(x, d, normalize=True)* 1./np.sqrt(n)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Before we begin, we should visualize our data matrix and verify that it satisfies the properties we want.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "n = 1000\n",
                "d = 63\n",
                "X = generate_X(n, d)\n",
                "\n",
                "# Show X and X^T X\n",
                "plt.imshow(X)\n",
                "plt.show()\n",
                "plt.imshow((X.T @ X))\n",
                "\n",
                "plt.show()\n",
                "print(\"diag(X^T X):\", np.diag(X.T @ X));\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "You should see that $X^T X$ is _approximately_ orthogonal. Think about why is it not exactly orthogonal. How would we change `generate_X` to ensure orthogonality?\n",
                "\n",
                "Now, we will write a function to generate the true weight vector $\\vec{w}$ that we are trying to learn, as well as the perturbed observation $\\vec{y} = X \\vec{w} + \\vec{\\varepsilon}$ that we can actually see.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "def generate_data(n, d, sigma=0.4, w=None, suppress_output=False):\n",
                "    X = generate_X(n, d)\n",
                "\n",
                "    # if w is not specified in the input, pick a\n",
                "    # random weight for each of the `d` dimensions\n",
                "    if w is None:\n",
                "        w = np.random.randn(d)\n",
                "        if not suppress_output:\n",
                "            print(\"true w:\", w)\n",
                "\n",
                "    # the standard deviation of the measurement error\n",
                "    eps = np.random.randn(n) * sigma\n",
                "    y = X @ w + eps\n",
                "\n",
                "    return X, w, y\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Now we've got some synthetic data. Let's look at the basic steps we need to take to optimize our hyperparameter:\n",
                " - Pick different values of the hyperparameter $\\lambda$\n",
                " - Train our model on some data using each particular value of $\\lambda$\n",
                " - Evaluate the performance of our model on *different* data\n",
                " - Pick the $\\lambda$ that yielded a model with the best performance.\n",
                "\n",
                "Notice one thing here - when evaluating the performance of our model, we evaluate it on *different* data that was not used in training. Why? Fundamentally, this is to prevent overfitting. When our model is trained on some input data, we can imagine that it tries to \"memorize\" the training data as best it can. If we then evaluated our model on the same training data, all we'd see is how well it memorized that data. In contrast, by evaluating it on *different* data, we are able to see how well it can generalize to new data it has not previously seen.\n",
                "\n",
                "Therefore, we will need to **split our dataset up into two parts: data used for training, and data used for validation.** You will do so here.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Generate a validation data set\n",
                "# Does it matter how we generate the split in this case?\n",
                "# How about in general?\n",
                "# Hint: what if the data source is sorted in some way?\n",
                "\n",
                "def generate_training_validation_split(X, y, validation_fraction = 0.2):\n",
                "    # Returns a tuple of 4 components:\n",
                "    #   X_train: (1-validation_fraction) rows of X\n",
                "    #   y_train: (1-validation_fraction) rows of y\n",
                "    #   X_val: the complement of X_train\n",
                "    #   y_val: the complement of y_train\n",
                "    n = X.shape[0]\n",
                "    n_validation = int(n * validation_fraction)\n",
                "\n",
                "    ### start Generate_Validation_Split ###\n",
                "\n",
                "    ### end Generate_Validation_Split ###\n",
                "\n",
                "    return X_train, y_train, X_val, y_val\n",
                "\n",
                "sigma = 0.6\n",
                "X, w, y = generate_data(1250, 63, sigma)\n",
                "X_train, y_train, X_val, y_val = generate_training_validation_split(X, y)\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "X_train.shape\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "def ridge_regress(X, y, lambd):\n",
                "    d = X.shape[1]\n",
                "    return np.linalg.solve(X.T @ X + lambd * np.eye(d), X.T @ y)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Let's first verify that our implementation of ridge regression is working. First, let's choose $\\lambda = 0$, so it reduces to basic least squares. **Use the training data to compute an estimate $\\hat{w}$ and plot it on the same axes as the true $\\vec{w}$**, to see if they are close.\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Show w vs w_hat\n",
                "### start Run_Least_Squares ###\n",
                "\n",
                "### end Run_Least_Squares ###\n",
                "plt.plot(w, label='w')\n",
                "plt.plot(w_hat, label='w_hat')\n",
                "plt.legend()\n",
                "plt.show();\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Now we can move on to hyperparameter tuning. One thing to notice is that we cannot actually choose $\\lambda$ to minimize $\\|\\vec{w} - \\hat{w}\\|^2$ directly, since $\\vec{w}$ is not known when working with real-world data. Instead, we should try to minimize $\\|y_{val} - X_{val} \\hat{w}\\|^2$ i.e. the error our model makes on the validation data. Intuitively, it is clear that these two quantities should be roughly proportional, and in our case when $X^T X = I$, it is easy to show that they are _exactly_ proportional.\n",
                "\n",
                "**Plot the validation error made after training the model using a range of lambdas. On the same graph, also plot the *true* estimation error $\\|\\vec{w} - \\hat{w}\\|^2$, on a separate vertical axis so the two plots are of similar scale.**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Show estimation error of w for a range of lambdas\n",
                "def plot_validation_errors(X_train, y_train, X_val, y_val, candidate_lambdas=np.logspace(-3, 0, 50), newplot=True):\n",
                "    w_mses = []\n",
                "    val_errors = []\n",
                "    true_errors = []\n",
                "    for l in candidate_lambdas:\n",
                "        ### start Compute_Errors ###\n",
                "\n",
                "        ### end Compute_Errors ###\n",
                "\n",
                "    if newplot:\n",
                "        plt.figure()\n",
                "        ax1 = plt.subplot(2,1,1)\n",
                "        plt.title(\"Error vs $\\lambda$\")\n",
                "        plt.yscale('log')\n",
                "        plt.xlabel(\"$\\lambda$\")\n",
                "\n",
                "        color = 'tab:red'\n",
                "        ax1.set_ylabel('Validation Error', color=color)  # we already handled the x-label with ax1\n",
                "        ax1.tick_params(axis='y', labelcolor=color)\n",
                "        ax1.plot(candidate_lambdas, val_errors, color=color)\n",
                "\n",
                "        ax2 = plt.twinx()\n",
                "        color = 'tab:blue'\n",
                "        ax2.set_ylabel('Estimation Error', color=color)  # we already handled the x-label with ax1\n",
                "        ax2.tick_params(axis='y', labelcolor=color)\n",
                "        plt.plot(candidate_lambdas, true_errors)\n",
                "\n",
                "        plt.xscale('log')  # Need to set scale after axes are created\n",
                "    else:\n",
                "        plt.plot(candidate_lambdas, val_errors)\n",
                "\n",
                "plot_validation_errors(X_train, y_train, X_val, y_val)\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "From the plot, you should see an obvious optimal value of $\\lambda$ that minimizes the validation error (if you don't, try regenerating the data matrix, you might have just been unlucky), and another value of $\\lambda$ that minimizes the true error. **Explain why these values are not the same**.\n",
                "\n",
                "Then, **write a function that numerically computes the $\\lambda$ that, when used to train a model on the training set, minimizes the validation error.***\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "_Your explanation here_:\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "def optimize_lambda(X_train, y_train, X_val, y_val, candidate_lambdas=np.logspace(-3, 0, 50)):\n",
                "    ### start Optimize_Lambda ###\n",
                "\n",
                "    ### end Optimize_Lambda ###\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "print(\"Best lambda = {}\".format(optimize_lambda(X_train, y_train, X_val, y_val)))\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Is this value the *exact optimal* choice of hyperparameter? No! Since we calculated it empirically from noisy data, it will itself be a noisy estimate of the optimal hyperparameter. It is important to note that *both* noise in the training data and, crucially, noise in the *validation data* will affect the \"optimal\" lambda that we obtain from validation.\n",
                "\n",
                "Let's see an example of this. In this next part, we will hold the training data fixed, and keep resampling our validation data from the same distribution. For each set of validation data, we will perform this hyperparameter optimization process again, and obtain a new \"optimal\" value of lambda.\n",
                "\n",
                "As a point of comparison, first **compute the optimal value of lambda.**\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "def compute_optimal_lambda(d, sigma, w):\n",
                "    ### start Optimal_Lambda ###\n",
                "\n",
                "    ### end Optimal_Lambda ###\n",
                "\n",
                "optimal_lambda = compute_optimal_lambda(d, sigma, w)\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "def repeatedly_optimize_lambda(n, d, validation_fraction = 0.2, num_samples = 50, sigma=0.4):\n",
                "    # first, we will compute a *fixed* set of training data, so the\n",
                "    # only source of randomness is from the validation set\n",
                "    X_train, w, y_train = generate_data(int(n * (1 - validation_fraction)), d, sigma, suppress_output=True)\n",
                "\n",
                "    plt.title(\"Error vs $\\lambda$\")\n",
                "    plt.yscale('log')\n",
                "    plt.xlabel(\"$\\lambda$\")\n",
                "    ax0 = plt.gca()\n",
                "    plt.twinx()\n",
                "    plt.yticks([], \"\")\n",
                "\n",
                "    optimized_lambdas = []\n",
                "\n",
                "    for i in range(num_samples):\n",
                "        X_val, _, y_val = generate_data(int(n * validation_fraction), d, sigma, w) # pass w in to keep it fixed\n",
                "        optimized_lambdas.append(optimize_lambda(X_train, y_train, X_val, y_val))\n",
                "        plot_validation_errors(X_train, y_train, X_val, y_val, newplot=False)\n",
                "        plt.twinx()\n",
                "        plt.yticks([], \"\")\n",
                "\n",
                "    # Have to set the scale last or the plots won't line up correctly\n",
                "    plt.xscale('log')\n",
                "    ax0.set_ylabel('Normalized Error')\n",
                "    return optimized_lambdas, w\n",
                "\n",
                "sigma=0.6\n",
                "optimized_lambdas, w = repeatedly_optimize_lambda(n=10000, d=d, validation_fraction=0.5, sigma=sigma)\n",
                "print(\"Mean optimized lambda: {}\".format(np.mean(optimized_lambdas)))\n",
                "print(\"Stdev of optimized lambdas: {}\".format(np.std(optimized_lambdas)))\n",
                "print(\"Optimal lambda: {}\".format(optimal_lambda))\n",
                "plt.show()\n",
                "\n",
                "plt.xlabel(\"$\\lambda$\")\n",
                "plt.ylabel(\"Frequency\")\n",
                "plt.hist(optimized_lambdas, bins=25)\n",
                "plt.show();\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Adjust `n` and `validation_fraction` to keep the size of the training set constant while varying the size of the validation set. The below code block might help with that.\n",
                "\n",
                "How do the accuracy and variance of the tuned lambdas vary as the size of the validation set increases? Why might the variance be large even with a large validation set?\n",
                "\n",
                "How do they change with the size of the observation noise, `sigma`?\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "n_train = 5000\n",
                "sigma=0.3\n",
                "for vf in [.1, .2, .4, .6, .8, .9]:\n",
                "    n_total = int(n_train / (1 - vf))\n",
                "    print('-' * 40)\n",
                "    print(\"Validation Fraction:\", vf)\n",
                "    optimized_lambdas, w = repeatedly_optimize_lambda(n_total, d, vf, sigma=sigma)\n",
                "    plt.show()\n",
                "    print(\"Mean optimized lambda: {}\".format(np.mean(optimized_lambdas)))\n",
                "    print(\"Stdev of optimized lambdas: {}\".format(np.std(optimized_lambdas)))\n",
                "    print(\"Optimal lambda: {}\".format(compute_optimal_lambda(d=d, sigma=sigma, w=w)))\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "_Your observations here:_\n",
                "\n"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": []
        }
    ],
    "metadata": {
        "kernelspec": {
            "display_name": "Python 3",
            "language": "python",
            "name": "python3"
        },
        "language_info": {
            "codemirror_mode": {
                "name": "ipython",
                "version": 3
            },
            "file_extension": ".py",
            "mimetype": "text/x-python",
            "name": "python",
            "nbconvert_exporter": "python",
            "pygments_lexer": "ipython3",
            "version": "3.8.5"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 4
}